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In this paper, we design a semi-implicit scheme for the scalar time fractional reaction- 
diffusion equation. We theoretically prove that the numerical scheme is stable without 
the restriction on the ratio of the time and space stepsizes, and numerically show that 
the convergent orders are 1 in time and 2 in space. As a concrete model, the subdiffusivc 
predator-prey system is discussed in detail. First, we prove that the analytical solution 
of the system is positive and bounded. Then we use the provided numerical scheme 
to solve the subdiffusive predator-prey system, and theoretically prove and numerically 
verify that the numerical scheme preserves the positivity and boundedness. 

Keywords: time fractional reaction-diffusion equation; subdiffusivc predator- prey system; 
positivity; boundedness. 
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1. Introduction 

Mathematically, the reaction-diffusion systems take the form of semi-linear 
parabolic partial differential equations. Usually, in real world applications, the re- 
action term describes the birth-death or reaction occurring inside the habitat or 
reactor. The diffusion term models the movement of many individuals in an envi- 
ronment or media. The individuals can be very small particles in physics, bacteria, 
molecules, or cells, or very large objects such as animals, plants. The diffusion is 
often described by a power law, {x 2 (t)) — (x(t)) 2 ~ Dt a , where D is the diffusion 
coefficient and t is the elapsed timePJ In a normal diffusion, a = 1. If a > 1, the 
particle undergoes supcrdiffusion, and it results from active cellular transport pro- 
cesses. If a < 1, the phenomenon is called sub-diffusion, it can be protein diffusion 
within cells, or diffusion through porous media. This paper concerns the subdiffu- 
sive reaction-diffusion system, which corresponds to the classical reaction-diffusion 
equation with its first order time derivative replaced by the a— th order fractional 
time derivative. 

As an important concrete example, we detailedly discuss the subdiffusive 
predator-prey model. All living things within an ecosystem are interdependent. A 
change in the size of one population or the environment they live affects all other or- 
ganisms within the ecosystem. This is shown particularly clearly by the relationship 
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between predator and prey populations. Cavani and FarkaP introduce diffusion to 
the Michaelis-Menten-Holling predator-prey model. Then the more general models 
are considerecP^IH. Here we further discuss the Michaelis-Menten-Holling predator- 
prey model with the subdiffusive mechanisirM 



d a N d N / AT aP \ 

— — = — -7T+N [l-N , xe(l.r), t>0, 

dt a dx 2 V P + NJ ' v ' h 

gap Q2p / 1 + 5 pp N \ 

erP -4; + 7^ — — ,ie(I,r), i > 0, 



(1.1) 



<9i a <9x 2 V 1+/3P P + N 



with the positive initial conditions and the homogeneous Dirichlet boundary con- 
ditions 

N(l, t) = N(r, t) = P(l, t) = P(r, t) = 0, (1.2) 
or the homogeneous Neumann boundary conditions 

(ON^X, £) /dx^ | x— l and r, respectively — [dP [x^ tjf 3xj \x=l and r, respectively — 0, (1-2 ) 

where a, er, and (3 are positive real numbers, and < 7 < 5. We prove that the 
analytical solution of (II. ip and ( (jl.2[) or (1.2')) is positive and bounded. 

For the analytical solution of the sub-diffusion equation, the reader can refer to 
Refs. [3 [H [9l [lOl [TTJ and the references therein. There are also some works for the 
numerical solutions of subdiffusion equations, e.g., Refs. [12l [13l HH [T5l [TBI HZ1 I n 
particular, Zhang and Sun develop the semi- implicit schemes for the subdiffusive 
equations^" 1 . For the past few decades, the semi-implicit schemes are widely used in 
various complicated time dependent non-linear equations. Usually the semi-implicit 
schemes use two time levels; in time level 1, the nonlinear terms are explicitly 
computed, and then to implicitly solve the high order linear terms. The expected 
advantage of the semi-implicit scheme is that as the nonlinear terms are computed 
efficiently but not losing good numerical stability. Here, we construct the semi- 
implicit scheme to numerically solve the subdiffusive reaction-diffusion equation. 
The stability of the numerical scheme is strictly proved, and it has no restriction on 
the ratio of the sizes of space steps and time ones. The convergent orders 1 in time 
and 2 in space are theoretically obtained and numerically verified. Moreover, we 
use the provided scheme to numerically solve the subdiffusive predator-prey model. 
We show both theoretically and numerically that it preserves the positivity and 
boundedness of the solutions of the subdiffusive predator-prey model. 

The outline of the paper is as follows. In Section 2, we propose the time fractional 
semi-implicit scheme for the subdiffusive reaction-diffusion equation. We discuss the 
stability and convergence of the proposed scheme in Section 3, and prove that the 
temporal approximation order is 1 and the order in space is 2. In Section 4, we first 
prove the positivity and boundedness of the solution of the subdiffusive predator- 
prey model, then certify that the numerical scheme preserves its positiveness and 
boundedness. In Section 5, we perform the numerical experiments to confirm the 
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convergent orders and positivity and boundedness preserving. We conclude the pa- 
per with some remarks in the last section. 

2. Scheme for the subdiffusive reaction-diffusion equation 

We first consider the following scalar subdiffusive reaction-diffusion equation: 

with x e fl = (0, 1), < t < T, < a < 1, the initial condition 

u(x,0)=g(x), xefl, (2.2) 
and the boundary conditions 

u(0,t) = u(l,t) = 0, 0<t<T, (2.3) 

or 

{du(x,t)/dx)\ x=0 = {du{x,t)/dx)\ x=1 = 0, < t < T, (2.3') 
where 9 g^f is the time fractional Caputo derivative defined as 
d a u(x,t) 1 /"* du(x,s) 1 



r du(x,s) 1 , n 

/ — tt^t r-ds, 0<a<l. (2.4) 

J ds (t-s) a V ; 



<9i Q T(l - a) J ds {t - s) a 

For ease of presentation, we uniformly divide the spacial domain ft = (0, 1) into M 
subintervals with stepsize h and the time domain (0, T) into N subintervals with 
steplength r. Let Xi — ih, i = 0, 1, • • • , M; tj — jr, j — 0, 1, • • • ,N. Let the grid 
function be u\ — u(xi, tj), < i < M, < j < N, denote 

and define 

D> * = r(2 - a) [< ~ - 6 ™->* - ( 2 - 6 ) 

as the discrete time fractional derivative^!, where bj = (j + l) 1- " — 2 l a - It can 
be noted that bj > and 1 = bo > &i > • • • > b n > (1 — a)(n + 1)~ Q . There exists 
the following error estimate between (d a u(xi,t)/dt a )\t=t n and D^uf: 

Lemma 2.1 (Ref. 18). Suppose < a < 1, and let u(xi,t) S C 2 [0,t n ], then 



d a u(xi, i) 



dt a 



n a n n 



<^T7T- 7- max \d 2 u(x u t)/dt 2 \ -T 2 ~ a . 

~ T(2 - a) o<t<t„ 1 v ' " 1 
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And it is well known that S^uj is the 2nd order central difference approximation of 
d 2 u(x,t)/dx 2 at (xi,tj). Replacing n by n + 1, Eq. (|2.6[) can also be recast as 

1— a j+1 _ u n-j 



rr-± — 

r ~r(2-a)^ 



1 n— 1 

— K +1 - - 6 j+1 )« w -' - M ]. 

1 — air" * — ; 



T(l - a). 



(2.7) 



Combining (|2.5[) and (|2.7|l . we design the semi- implicit finite difference scheme of 
(|2~T]) as 

I?«C/f +1 = 8 2 X U? +1 + f(U?), (2.8) 

where i = 1, 2, • • • , M — 1 for boundary condition (|2.11[) . and i = 0, 1, • • • , M for 
(2.13'), n — 1, 2, • • • , N — 1. Denoting T(2 — a)r a by C Q we rewrite the above 
semi-implicit finite difference scheme as 

E b i~ = C « 5 l U ? +l + C a f(U?). (2.9) 

3=0 

From (|2.2I) the initial condition is specified as 

0? =$(-*), for i = 0,l,.-. ,M; (2.10) 
and from (|2.3I) or (2.3') the boundary conditions are given as 

E/? = C/& = 0, for n = 0,l,--- ,JV, (2.11) 

or 

U\ = U{\ U n M+1 = U^_ 1} for n = 0, 1, • ■ • ,N, (2.13') 
namely, the central difference discretization is used for the Neumann boundary. 



3. Stability and convergence of the numerical scheme fl2.9p ~ fl2.llf> 

Now we discuss the stability and convergence of the numerical schemes, first we 
analyze the numerical stability. Let Z7™ be the approximate solution of the numerical 
scheme (|2~9l - (|2~TTl ). and denote e™ = J7f - UJ 1 . From (T231) . we immediately obtain 

n n 

= E ' - E b ^ +1 ~ J + ( £ m - k +1 +^)+c a (m)-m)), 

(3.1) 

while the perturbation errors of boundary conditions are 

eo = 4 = 0, 1 < n < AT. (3.2) 
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Throughout the paper, we assume that the function / satisfies the local Lipschitz 
condition with the Lipschitz constant L, namely, 

l/OO,*)) - f{u{x,t))\ < L\u(x,t)-u(x,t)\, 

when \u{x,t) — u(x,t)\ < e$, x £ il and t 6 [0,T], where £o is a given positive 
constant. 

Denote e" = (e™,e2, ■•• ,^i_i) T - We define the discrete L 2 norm as ||e"|| = 

M \V2 

^ S e ? ) ■ Numerical stability result is as follows. 
Theorem 3.1. The numerical scheme h2.9fy - l2.11\) is stable and there exists 

when T < (l/(r(l - ajL)) 1 /". 

Proof. We prove this theorem by mathematical induction. Taking n = in (|3.ip . 
we have 

(* + ^r) 4 = + + ej_ x ) + C Q (/([/«) - f(Ui))- (3.4) 

Multiplying both sides of (I3.4[) by he] and summing up, there exists 



1 1 1 2 



Af-1 „ M-l M-l 

= £ (e? • 4 ■ h) + ( (e i +1 + eU) ■ el ■ h) + C a £ {f(Uf) f(U?))e}h 

i—l i—1 i—1 

or 1 

<l|e ||-|| e 1 || + ^||e 1 || 2 + C Q .L.||e°||.||e 1 ||. 
Then we obtain 

He 1 !! < (l + <7 a £)[[e°||, 

it can be easily checked that this means (|3.3[) holds for e l . Now supposing (|3 . 3[) 
holds for e , e 2 , • • • , e™, we prove 

" e " +1 ^ (l-a)-^r(2- a )L " e °"- 
Just as the above process, multiplying both sides of (|3.ip by /ie" +1 and summing 
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up, we get 

-i- . 

h? 

n M-l n M-l 

= E E w j <r l h -EE b^-'^h 

j—Oi—l j—1 i—1 

„ M-l M-l 

i=l i=l 
ra-lM-1 M-l 

= E E & - Wer'e?^ + E bn^h 

j=0 i—1 i—1 

~ M-l M-l 

i=l i=l 

< $> - ^+i)ll« n - i ll ■ l! e " +1 H + M^]] ■ l! en+1 H + ^He n+1 || 2 

+ C Q -L.||e"||.||e n+1 ||. 
Then, there exists 

< ||e°|| 

" (1 - a) -T«r(2-a)L" 11 



□ 



Theorem 3.2. Let u{xi,t n ) and XJf be the exact solutions of the subdiffusive 
reaction- diffusion equation J,g.j) )- fO)) and of the numerical scheme I12.9\ )- ^2JJ]) . 

respectively, define e" = u(xi,t n ) — U\ l and E n — (e™,^ - ' - ) £ m-i) T ? ^ erl 
E n satisfies the following error estimate: 

\\E n \\<C(r + h 2 ), (3.5) 

when T < (1/(T(1 - a)!,)) 1 ' 01 . 



Proof. According to equation (|2.8[) . we know 

D?U? +1 = 8lU? +1 + f{U?). (3.6) 
By Lemma 2.1, there exists a positive constant C, such that 

< Ct 2 "", (3.7) 



d a u(Xi,t) n a„,Ti+l 
7^ | t=tn+1 -^ r M, 



I 



I 
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and it is well known that 



d 2 u(xi,t 



n+l) r2„ n+1 



dx 2 



Mi 



< Ch 2 . 



Thanks to the Lipschitz continuity of / with respect to u, we have 



(3.8) 



<L < +i -< <Ct. 



(3.9) 



Based on (I3.7I)-(I3.9I). there exists 



where Rf — 0(t + h ). For convenience, denoting C Q i?" as r™ , we get 



V n J 3=0 j=l 

+C a (f(u(x i ,t n ))-f(U?))+r?+ 1 



(3.10) 



Similar to the proof of Theorem 13. 1[ taking n = in p.lOj) . multiplying both sides 
of (|3.10j) by he}, summing up and directly calculating, there exists 

H^H 2 < b \\E°\\ ■ H^ll + C a L\\E°\\ ■ 1 1^1 1 + Hr 1 !! • H^H, 

where r 1 = (r\, ■ ■ ■ ,t\f_- i ). Then 



I^H < (l + CaLJH^II + Hr 1 ! 



< 



bo 1 



(l-a)-T a T(2-a)L 



(\\E°\\ + ||r 1 ||). 



Now suppose that 



P m ||< 



b m -i 



(1 - a) - T a T(2 - a)L 



(\\E° 



l<j<m 



II). 



(3.11) 
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holds for m = 1, 2, • ■ • , n, then we prove that it holds for E n+1 . Similarly, 



1 + -^ \\E n+1 \\ 2 



2Ca 
h 2 

M—l n M—l n 

E E^r^- E E 6 i £ " +W£ ? +1/i 

i—l j—0 i—1 j=l 

M-l ^ M-l 

+ E -J^( s iti+ £ i-i)^ +1 h + C a E {f{u{xM)-f(U?))e r l +1 h 

i=l i=l 



E r? +1 e? +1 h 



n-1 



+ C a L\\E n \\ ■ \\E n+1 \\ + \\r n+1 \\ ■ \\E n+1 \\. 
Dividing by ||i? n+1 || at both sides, we conclude that 

n-1 

< £(&; - + M£°H + C a L\W n \\ + ||r" +1 ||. 

3=0 

According to p. lip , combining with bj 1 < bj +1 , we have 
\\E n+1 \\ 

^ g (i- Q )-^r(2- Q) L (l|i;; 11 + i^-i l|r ll} 
S (i-o)-r.'r(2- a .)L (l|E ° ll + ^ ll ' J|l) - 

Notice that 

£J? = 0, J5g = £^ = 0, < i < M, 1 < j < n + 1. 
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Together with b' 1 < and (n + l)r < T, we obtain 

ll^" +1 H<7 x v —, rr max Mr 3 !! 

11 11 ~ (1 - a) - T a T(2 - a)L i<j<n+i 11 11 

< m ax ILR'N 

~ (1 - a) - T a T(2 -a)L i<j<n+i 11 11 

T Q r(l - a) 
" (l-a)-T°T(2-a)L i<?<n+i 

<C(r + ^ 2 ). □ 

Remark 3.1. The above analysis focuses on scalar equation, but it can be easily 
extended to the vector one. Take Ijl.ljl as an example, and denote the nonlinear 
term of the first equation as f(N, P) and the nonlinear term of the second equation 
as g(N, P). Denoting p n = N(x z , t n ) - N™ and rf = P(x l ,t n ) - P z n , and using the 
following simple trick 

\f(N? +1 ,P? +1 )-f(N?,PP)\ 

= i/(iv-+\ff +i ) - f{N? + \p?) + f(N?+\pn - f(N?,pn\ 

< |/(iVf+ 1 , J P^+ 1 ) - f(N?+\Pn\ + \f(N? +1 ,PD - f(N?,P?)\ 
<L(\P? +1 -P?\ + \N? +1 -N?\), 

lead to 

II^H ^ 7i — r^kko — u (^°» + ll^ " + ^ H^D- 

(1 — a) — 1 a i (2 — a)L i<j<»™ 

and 

I' 7 '" 1 " ^ 7i ^ W79 T7 ^ + m + ^ a < x l|rJ|l) ' 

(1 — a) — 1 a L (2 — a)L i<j<m 
for m = 1, 2, • • • , n, by the analysis similar to Theorem 3.2. 

Remark 3.2. The above analysis is for L 2 estimates, and we hope that the H 1 
estimates can be obtained in the near future, but it needs more delicate tricks to 
dealing with the no nlin ear terms. In fact, there are already the H 1 estimates for 
the linear equations 

4. Positivity and boundedness of the analytical and numerical 
solutions of the subdiffusive predator-prey model 

In this section, we first prove that the analytical solutions of (|l.l[) - (ll.2j) are positive 
and bounded, then demonstrate that both the numerical schemes (|2.9p - (|2.10p with 



15, 2013 2:24 



Revised'SCM-2012-0423 



10 Yanyan Yu, Weihua Deng, Yujiang Wu 

(|2.11l) and with (2.13') preserve their positivity and boundedness when utilized to 
numerically solve (|1. lft - (11 .21) and (II . If) and (1.2'), respectively. 

4.1. Maximum principle for analytical solutions 

For analyzing the properties of the analytical solutions, we introduce the following 
maximum principle. The considered equation is 



Lu = — - 7 — + c(x,t)u = f(x,t), OM)6fi T = ftx(0,T], (4.1) 



where £1t is a bounded domain with Lipschitz continuous boundary. 

Theorem 4.1. Assuming that the coefficient c(x,t) > and f(x,t) < (resp. 
f(x,t) > 0) in Q,t and u € C 2,1 (r2r) p| C(f2y) is the solution of J^. lj) , then the 
non-negative maximum (resp. non-positive minimum) value of u(x,t) in Qt (if 
exists) must reach at the parabolic boundary Ft, namely 

maxu(i, t) < m&x{u(x, f), 0} (resp. wmu{x 1 t) > min{u(x, t), 0}). (4-2) 

In fact, if the non- negative maximum value of u(x,t) is not at the boundary Tt 
and f(x,t) < 0, then there exists a point (x*,t*) € fix such that 



u{x* , t*) > max{w(x, t), 0} and u(x*,t*) > maxu(x,t). 
Ft n T 



Let b > 0, for any e > 0, we introduce the auxiliary function 



v(x, t) 



u(x 7 t) - ee . 



On the one hand, we know that, for any {x,t) G fix, v(x,t) satisfies 




(4.3) 



f(x, t) - e(&i 1 -« J B 1)2 _ a (W) + e bt c(x, t)) < 0, 



where E a ^(z) is the Mittag-Leffler function. At the maximum point (x*,t*), ac- 



15, 2013 2:24 



Revised'SCM-2012-0423 



Positivity and boundedness preserving schemes for the fractional reaction- diffusion equation 11 

cording to the definition of Caputo derivative, we have 
d a u(x*,t) 1 f du{x*,s) 1 



dt a r(l - a) J Q ds (t* - s) a 

= l[m n r(o - ; Y,bj-Wx*,t*-jT) - u(x*,t* - (j + l)r)), 
t->0 LIZ — a) *■ — ' t 

= lim n r/o " - C 1 - bi)u(x*,t* - r) b n u(x' , 0)) 

= nn ), m T^ 1 -&i)( u O ,* ) - -r)) + •■• 

t->0 L (Z — a) 

+ b n {u{x*,t*)-u(x*,0))) 

^ lim n rfl ~ s b n (u(x*,t*)-u(x*,0)) 
r-s-0 1 (2 — a) 

>hm — (u(x ,t)-u(x ,0 

r->0 1 (2 — a) 

(1 - a)T~ a 
> { r(0 ' , (u(x*,t*)-u(x*,0)), 
1 (2 — a) 

where r = t*/n, bj is defined in ()2.6)) . and 1 = 6o > b\ > ■ ■ ■ > b n > (1 — a)(n + l)~ a 
is used. Since u(x*,t*) > u(x*,0), denoting m* = u(x*,t*) — u(x*,0), there exists 



-ds 



dt a 1 dt a 

when e < r ( 2-afel-^"- a (^-) - Together with = 1 ' 1 ^ 

know 

which is contradictory with (|4.3p . Similar analysis can be done for the case f(x, t) > 
0. So, from the above analysis, we arrive at Theorem 14. II 

Remark 4.1. If we take the initial condition of (14.11) as u{x, 0) = 0, and its bound- 
ary conditions are Dirichlet's and homogeneous, then from Theorem 14.11 we have 
u(x,t) < when f{x,t) < and u(x,t) > when f(x,t) > 0. The same results 
still hold if the homogeneous Neummann boundary conditions are used, since the 
maximum (minimum) value of u(x, t) at the boundary is non-positive (nonnegative) 
under the homogeneous Neummann boundary conditions. In fact, if the maximum 
value of u(x, t) at the boundary is positive, suppose it is located at the left boundary 
(similar analysis can be done if at the right boundary) and denote the one closest 
to the line t — by u(xi,t*), then for any given sufficiently small e, there exists 
i E e (xi,xi + e) such that u(x t + e,t*) - u{x u t*) = u ^f ] \ x=ic )e 2 < since 
^"fof \x=xi — 0. So there exists sufficiently big M > 2 and small St > such that 
< for any (x, t) e fi* = {(x,t) \x t < <x< Cs+jj < x t +e; t*-St < 
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t <t*}. Now consider (14.1j) in the domain il* = {(x, t) |f e < x < xi +e; < t < t*}, 
obviously the maximum value of u(x, t) still is obtained at the parabolic boundary 
r t * of tt* . Furthermore, if taking e small enough, then the maximum value is lo- 
cated in the domain r t » n Q*. Now at the maximum point, O^g- — ^4? + cu > 0, a 
contradiction is reached. 



4.2. Positiveness and boundedness of the analytical solutions 

Using the upper and lower solutions method, we prove the positiveness and bound- 
edness of the analytical solutions (ll.l[) - (|1.2p and (jl.ip and (1-2'). First, we introduce 
the definition of the upper and lower solutions. 

Definition 4.1. For the system of equations (i — 1, 2) 

(Oil ' \ 
Uior-^j = gi (x,t), x€dtt, t€(0,T], (4.4) 

m(x, 0) = ifi{x), x e fi, 

suppose that Ui(x,t) and Ui(x,t) satisfy 

Bui-gi(x,t) > > Bm -gi(x,t), x G 90, t £ (0,T], (4.5) 
Ui(x,0) — ifii(x) > > Ui(x,0) — ifi(x), igSI, (4-6) 

and /i(-, ■) is quasi-monotone decreasing, f 2 {-, •) is quasi-monotone increasing, and 
ga ~ q2~ d a ui d 2 m 

■W ~ J? ~ * ° >-^r - - A(«..«). («) 



2,-, 9 Q u 2 d 2 u 2 



f2(u u u 2 )>0>—^ ^-f 2 (ui,u 2 ), (4.8) 



dt a dx 2 ' 9i™ 9a; 2 

then U(x,t) = {u\{x,t),u 2 {x,t)) and V(x,t) = (u\{x,t),u 2 {x,t)) are respectively 
called upper solution and lower solution of the system (|4.4j) . 



Theorem 4.2. Suppose {fi,f 2 } is mixed quasi-monotonous and Lipschitz contin- 
uous with respect to u\ and u 2 

\fi(ux,u 2 ) - fi(vi,v 2 )\ < L(\ut - Vx\ + \u 2 - v 2 \), 

where L is constant. If the upper and lower solutions, U{x,t) and V(x,t) 7 satisfy 
V(x,t) < U{x,t), then CT^P has a unique solution in \V{x, t), U(x, t)] . 



Proof. See Appendix A. □ 

Let us denote f(N,P) = N(l- N- -gfa) and g(N, P) = a P{-^f- + p^). 
We certify that both (|l.ip - (|1.2j) and (11.11) and (1.2') have lower and upper solutions, 
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(0,0) and (l,£i), where L\ > 1/7. First, it is needed to check that / is quasi- 
monotone decreasing function and g is quasi-monotone increasing function. Since 
the derivatives of the nonlinear terms are 

df(N,P) _ aN 2 
dP ~ ~(P + N) 2 

and 

dg{N,P) _ aP 2 
ON ~~ (P + N) 2 ' 

it is obvious that / is decreasing w.r.t P and g is increasing w.r.t N. At the same 
time, it can be noted that both (|1.2[) and (1.2') satisfy the conditions (|4.5p directly. 
And it can be easily checked that the upper solution U(x,t) = (N(x,t), P(x,t)) = 
(l,£i) and the lower solution V(x,t) — (N(x,t),P(x,tj) = (0,0) satisfy 

d«N d 2 N - d ° N 92 N 

gap Q2p d a P d 2 P 



6* dx 2 ^^^--^-gUP). (4.10) 

So if we specify the initial condition of (|l.ll) such that N(x, 0) € [0, 1] and P{x, 0) € 
[0, Li] for any x £ (l,r), then the initial condition satisfy (|4.6[) . From Theorem 14.21 
the exact solution of (|1.1[) - (|1.2I) (or (jl.ip and (1.2')) is bounded and positive. 

4.3. Positiveness and boundedness of the numerical solutions 

We show that the numerical schemes ([279])- (2. 13') preserve the positiveness and 
boundedness of the corresponding analytical solutions of (|l.ll) -f 1.2'). 

First, for (jTTTj) with the initial conditions N(x,0) £ (0, 1] and P(x,0) e (0,ii], 
L\ < ~, for any x £ (lff)> we have its discretization scheme 

n-1 

- c a 8tNr x = - f>j+i)Nr j + 

j'=o 

+c Q ivr(i-7vr 



pn + N n, 



n-1 

p? +i - c a s 2 x p? +i = Y,(bj - b j+ i)pr j + b n p? 
3=0 



ry + §fjpn N n 

n a p n ( 1 J— 1 4. ill ) 
a 1 1 1 + pPr 1 P? + N?'' 



We use the induction method to prove < NJ 1 < 1 and < P" < L\ for any i and 
n. First, < N° < 1 and < PP < L\ hold obviously. Now suppose < Nf < 1 
and 0<P l k < Li for any k < n, we prove that it still holds when fc = n + 1. 
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First of all, denote w :— C a N™(l — N™ — p °^„ ) and v := C a aPp(~^ff- 
pn+N n ) • When C a < 1, it's easy to obtain 



and 



w < C„JV?(1 - Nl l ) <C a (l- N?), (4.11) 
w>C a N?(l-N? -a) > -C a aN?; (4.12) 



V < C^m—^ + ^n)< Ca*(l - ^7), (4.13) 

v > C a aP?{-5 + pn N J Nn ) > -C a aP?8. (4.14) 



Because of bj > b j+ i, < TV™ < 1 and < Ff < L x , we know 



(i - w < - ^+i)^r _j + &»^° < + (4-i5) 



3=0 

and 

n-i 



(1 - h)P? < J2(b 3 - b j+ i)Pr j + bnP? < + ( 4 - 16 ) 

3=0 

Owing to (tlTTl) . dHHl) and P~T5]1 . when C a < min{±^, ±=£±}, we get 

< N? +1 - C a S 2 x Nl l+1 < 1. (4.17) 

In the similar way, owing to (|4.13[) . (|4. 14[) and (|4. 16[) . when C a < — J'^ , we get 

< if + 1 - C a %P? +1 < Li. (4.18) 

From (|4TT7|) and (|4TT5|) , we can get < iV™ +1 < 1 and < Jf +1 < L x . In fact, if 
< N" +1 < 1 doesn't hold, then there exists i such that 

7V™ +1 < or iVf +1 > 1. 

If N" +1 < 0, then we choose the minimum in i = 0, • • ■ , M, and denote it by N£ +1 , 
which is non-positive. Thanks to (|4.17l) . we know that 

A^ +1 " ^(N^l - 27V«+i + iV^ 1 ) > 0. 
So 7V™ +1 < implies 

jyn+1 ^ fc— 1 fc+1 

fe 2 

Then we get N™+1 < N r k l+1 or N^+l < N™ +1 , which still hold even at the boundary, 
including the Dirichlet and Neummann boundaries. This is contradictory with the 
assumption that N£ +1 is the minimum. If N" +1 > 1, then choose the maximum 
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in i = 0,-- - ,M, and denote it by iV™ +1 , and N? +1 > 1 holds. Since iV™ +1 - 
%-{N?+i - 2iV" +1 + < 1, then 

^-(N l n + 1 - 2iv ; n+1 + Nf+i) > n; i+1 - 1 > 0. 

So 

. h AT,™ 4 : 1 + ATTVV 1 

jyn+1 - t — 1 1 + 1 

' 2 

Then we get A^™ +1 < A^ 1 or iV" +1 < 2V^\ which still hold even at the boundary, 
including the Dirichlet and Neummann boundaries. This is contradictory with the 
assumption. 

Similarly, we can verify that < P[ l+1 < L\. 

5. Numerical experiments 

We present the simulation results of the schemes (|2.9j) - (|2.11j) for Dirichlet bound- 
ary and (|2.9[) - (|2.10[) and (2.13') for Neummann boundary to verify all the above 
theoretical results. In particular, the subdiffusive predator-prey model (|1.1[) with 
homogeneous Neummann boundary conditions (1.2') is simulated, and the pictures 
are displayed. Example 5.1 and 5.2 numerically confirm the unconditional stability 
of the numerical schemes and first order convergence in time for any a € (0,1). 
Example 5.3 is for the subdiffusive predator-prey model with specified initial and 
boundary conditions. 

In the computations of Examples 5.1 and 5.2, we take the spacial steplength 
h = 0.0005, which is small enough so that the spacial error can be neglected for 
obtaining convergent rate in time direction. The errors are measured at time T = 1 
and by l°° norm. And a is, respectively, taken as 0.3, 0.6 and 0.9. 



Example 5.1. For (|2.1j) - (|2.3j) . we take its exact analytical solution as 

u(x,t) = t 2 siii(2Trx), (5.1) 

and the non-linear term as 

f(u) - — !— . (5.2) 
Then, on the right hand side, we need to add the forcing term 

,(*,*) = j^^sin M +4*^(2**) - t2s . m( l x)+4 , (5-3) 
and the corresponding initial and boundary conditions are respectively 

it(x,0) = 0, (5.4) 



u(Q,t) = u(l,t) = 0. 



(5.5) 
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Table 1 . The error and convergent rate of the 
proposed scheme for Example 5.1, when a = 
0.3 and h = 0.0005. 



e(h,r) 



128 
256 



1.628403729893591e - 003 

8.461705972403477e - 004 0.94444 

4.303180907271331e - 004 0.97555 

2.166950795682299e - 004 0.98974 

1.088150760031326e - 004 0.99379 

5.474834928109740e - 005 0.99099 



Table 2. The error and convergent rate of 
the proposed scheme for Example 5.1, when 
a = 0.6 and h = 0.0005. 



128 
256 



e(h,r) rate 

2.021632915774174e - 003 

9.783209560871864e - 004 1.0471 

4.709883249863767e - 004 1.0546 

2.270387164748922e - 004 1.0528 

1.100565972464995e - 004 1.0447 

5.382548524024422e - 005 1.0319 



Table 3. The error and convergent rate of 
the proposed scheme for Example 5.1, when 
a = 0.9 and h = 0.0005. 



i 

128 



e(h,r) rate 

3.392154303510031e - 003 

1.657079757972468e - 003 1.0336 

8.035388508833563e - 004 1.0442 

3.883816156070585e - 004 1.0489 

1.876360100631080e - 004 1.0495 

9.083194250991689e - 005 1.0467 



Example 5.2. For (|2.1|) . (|2.2[) and (2.3'), the exact solution and the boundary 
condition are, respectively, taken as t 2 cos(27rx) and 

du(x,t) _ du(x,t) _ 

We still use (|5.2j) and (|5.4|) as the non-linear term and initial condition, respectively. 
And the following forcing term is needed to add to the right hand side of the 
equation, 

9{x, t) = 2 / 2 ~ a) cos^Tra) + AirH 2 cos(2^) - } (5.6) 

1 (3 — a) t z cos(27ra;) + 4 



i 



i 
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Table 4. The error and convergent rate of the 
proposed scheme for Example 5.2, when a = 
0.3 and h = 0.0005. 



128 
250 



e(h,r) 



2.887575706533641e 


- 003 




1.527198871648983e 


- 003 


0.91897 


7.855525404816266e 


- 004 


0.95911 


3.983144175851994e 


- 004 


0.97980 


2.006793113920047e 


- 004 


0.98902 


1.009537078571210e 


-004 


0.99120 



Table 5. The error and convergent rate of the 
proposed scheme for Example 5.2, when a = 
0.6 and h = 0.0005. 



r e(h,r) rate 



2.713892842113319e 


- 003 




1.337146145350632e 


- 003 


1.0212 


6.550178686808295e 


- 004 


1.0296 


3.205477251759792e 


- 004 


1.0310 


1.572572233041747e 


- 004 


1.0274 


7.755217692539951e 


- 005 


1.0199 



Table 6. The error and convergent rate of 
the proposed scheme for Example 5.2, when 
a = 0.9 and h = 0.0005. 



t e(h, r) rate 



1 


3.769255105664726e 


- 003 




1 


1.832929458344568e 


- 003 


1.0401 


8.885772877893494e 


- 004 


1.0446 


128 
256 


4.302894303225280e 


- 004 


1.0462 


2.084713949954686e 


- 004 


1.0455 


1.012302249301378e 


- 004 


1.0422 



Example 5.3. Consider the reaction diffusion equation^ 

d 2 N n P 



d 2 P 7 + SpP N 

8x 2 [ 1 + 0P P + N 

with the homogeneous Neummann boundary conditions on the domain Q = [0, 1] 



Let us denote f(N,P) = N(l - N - g(N,P) = a P{-^§f + p^), 

and define (N, P) as the equilibrium point of (|5.7[) . In the case, a — 1, a = 1.1,7 = 
0.05, ft = 1, and S = 0.5, then as f(N,P) = and g(N,P) = 0, we can obtain the 



i 



i 
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equilibrium point (N,P) = (0.113585,0.471397). The simulations were performed 
for the system on a fixed grid with spatial stepsize h — 0.005 and time stepsize 
r = 0.1. As the initial condition, we use 

N(x, 0) = N + 0.0214 cos(tt:e), 
P(x, 0) = P + 0.0066 cos(ttx). 

We focus predominantly on displaying the properties of the numerical solutions 
for different time fractional order a, see Fig. 1, Fig. 2, and Fig. 3. 

6. Conclusions 

We introduce the unconditional stable semi-implicit numerical schemes for subdiffu- 
sive reaction diffusion equation with Dirichlet boundary condition and Neummann 
boundary condition, respectively. And the subdiffusive predator-prey model is de- 
tailedly discussed. We prove that its analytical solution is positive and bounded. 
Then we show that the proposed numerical schemes preserve the positivity and 
boundedness of the analytical solutions. The extensive numerical experiments are 




(a) N (b) P 



Fig. 1. Numerical solution for di = 0.005, d,2 = 0.2, a = 0.2. 




(a) N (b) P 



Fig. 2. Numerical solution for di = 0.005, d 2 = 0.2, a = 0.5. 
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performed to confirm the theoretical results and show the dissipative properties of 
subdiffusive predator-prey model. 
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Appendix A. 

Proof of Theorems 4.2: Taking the initial iteration function as 

,_(0) _(0)n /- - \ 

i (°) (°h t \ 
KM ,M ) = (^i'^2j, 

with u\ ^ Hi and u 2 ^ U2, define the following iteration 



d a u (k) 


d 2 u (k) 


dt a 


dx 2 


d a u {k) 


d 2 u (k) 


dt a 


dx 2 


S a uf 


d 2 u[ k) 


dt a 


dx 2 


d a ^ k) 


d 2 u 2 k) 




dx 2 


Buf ] \ da> 


;(0,T] = 




= m ( 



J _(fc) J -(fc-1) 



(fc-1) Jk-l), 



L ■ u 



(k) 



(fe-i) 



L-u 



(fc) 



L ■ u 



f2(u\ 



(fc-1) -(fc-1), 



/i(ui' 



(fc-1) -(fc-1) 



i "2 



(fc-i) , (fc-i) 



_ n„(*)i 



(A.l) 




(a)N 



(b)P 



Fig. 3. Numerical solution for di = 0.005, d 2 = 0.2, a = 0.9. 
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where L is the maximum of Lipschiz constants of /i and fi- Subtracting (14.7[) and 
(|4~8f from (|A.1|) leads to 



^ + L-(^ 1) -4 0) )<0 







a 2 (w^ - 




dt a 




dx 2 




d a (u ( 2 1] - 


4 0) ) 


d 2 {u 2 1] - 


4 0) ) 


dt a 




5a; 2 




d a (u^ - 




a 2 (ul 0) - 




dt a 




dx 2 




d a (u 2 0) - 




d 2 ^ - 




dt a 




dx 2 





+ ^-(nl° ) -ul 1) )<0 

+ L.(u^-u«)<0 

B(uj x) -u! 0) )|anx(o,T] = B(uj 0) -u?%nx { o,T] = 0, (i = 1,2) 
k uj 1) (a:,0)-tij 0) (a:,0) = u| 0) (x, 0) - u^ (x, 0) = 0, (a 6fl,i= 1,2). 

According to the maximum principle Theorem 14.11 and Remark I4.1[ we know that 



Supposing u\ < u[ k , u 2 - fe ^ < u^, and then that the nonlinear term fx is 
quasi-monotone decreasing results in 

h{ui ,u 2 J - Jn u i >U2 ) 



Ji( u i 7H2 ) - Jn u i 7H2 ) + h( u i (Ha J - /H u i 



(fc) ffc-lh 



In a similar way we also get 

J2{U 1 ,U 2 ) — J2\U 1 ,U 2 j < L ■ (U 2 



-(fe-1) „(fc-l) 



(A.2) 



(A.3) 



Jl Uil ' 2 / X 1 \— 1 ; u 2 ) S *-> ■ Oil 



(fc-1) -(fc-lh 



(fe) _ (fc-1) 



(A.4) 



(A.5) 
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So, together with (|A.2[) - (|A.5|) . the iteration (|A.1|) implies 

aaf.-Xk+l) (k), o 2 /-(fc+l) -(fc)s 

< iv • [u x ) + M% ,u 2 ,u 2 J < 1) 

Q Q /-(fc+l) -(fc)% ^^(fc+i) 

a (u 2 - u 2 ) a (u 2 - u 2 ) ,-(k+i) _(fe)x 

+L ' {U2 ~ U2 > 

< L ■ (4 fe) - ut 1] ) + h(4 k \^) ~ Huf^Mt 1 ") < 
yfczM-uT 1 )) ^(uf»-uf +1) ) (fc) cm.!, 

< L ■ (u[ k ^ -uf») + ACuM ) - ACui*- 13 .^) < o 

aof,.W n ,( fc+1 )A ^/..(fe) „(*+!) <\ 

+ L ' ( - 2 ~- 2 } 

S 1j ■ iu 2 - u.2 J + J2 Uii ,u 2 j - /aUii ) S U 

- uf ) )| 9 nx ( o,T] - B(u, (fc) - u! fc+1) )|a^(o,T] = 0, (i = 1,2) 

_ 0) - u< fc) (;r, 0) = uf } (or, 0) - uf 0) = 0, (z G ft, t = 1, 2). 

Then there exists 

< zf\ uf < uf +1) (fe = 0,l,.-0- 
Recalling the iteration (jA.lj) again, we deduce 

< l • (h£- 1 > - s?- 1 )) + /xCuf- 1 ),^- 1 )) - A^r 1 ),^- 1 )) < o 

g a (4 fc) -4 fc) ) ^(u£-uf) w w 

^ 9^ +L ' [ ~ 2 ~" 2 } (A.6) 

< A, • (U 2 - U 2 ) + >u 2 ) - J2(U 1 ,u 2 ) S U 

B(ul fc) -^' s) )laax(o ) r]=0, (t = l,2) 
I u, (fe) (a;, 0) - ui fc) (a!, 0) = 0, (x G fi, t = 1, 2), 
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then 



So 



Ui < u W < 



(fe) , _(fc) 



, (fe) , _(fe) , , _(i) 

< U. - < U - < ■ ■ ■ < U - 



<u u (i = l,2). 



Note that /i is quasi-monotone decreasing and /2 is quasi-monotone increasing, 
then there arc 



lim u\ 

fe^+oo 



(fe) 



Ui{x,t), 



lim uf } =^(3,*), (i = 1,2), 



which satisfy m > and 



' d a U! 


d 2 u\ 


dt a 


dx 2 


d a u 2 


d 2 u 2 


dt a 


dx 2 




d\i 


dt a 


dx 2 


d a u 2 


d 2 u 2 


dt a 


dx 2 



/i(ui,u 2 ) = 
/ 2 (wi,u 2 ) = 
/i(ui,w 2 ) = 

/2(Ul,U 2 ) = 

Bu^onxio.T] = £uJanx(o,T] = 
_ Ui(x,0) = Ui(x, 0) = </?i(x), (a; e fi, i = 1,2). 



Next we will certify 



1,2). 



Defining wi = u\ — u 1; w 2 — u 2 — u 2 , we nave known wi > and w 2 > from 
above discussions. According to the iteration, we can obtain 



d a u>i d 2 W\ 



dt° 



dx 2 



= /i(wi,u 2 ) - /i(ui,u 2 ) 

= /l(Ml,U 2 ) - /l(U!,U 2 ) + /l(Ui,U 2 ) - /l(U!,M 2 ) 

< L ■ (ui — u 1 ) + L ■ (v,2 — u 2 ) 
= L ■ (wi + w 2 ), 



d a W 2 3 2 W2 



dt° 



dx 2 



/ 2 (U1,M 2 ) - / 2 (U!,U 2 ) 
= /2(U1,U 2 ) - / 2 (U1,U 2 ) + / 2 (Ml,U 2 ) - / 2 (Ui,U 2 ) 

< L ■ (u 2 — u 2 ) + L ■ (ui — u^) 
= L ■ (u>i + w 2 ). 
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So 



r d a w 1 


d 2 wi 


dt a 


dx 2 


d a w 2 


d 2 w 2 


\ dt a 


dx 2 



- L ■ (wi + wi) < 



Bwi = U,tfw 2 = U 
Wt(x, 0) = 0, w 2 (x, 0) = 



0. 



Denoting w — W\ + w 2 > 0, it's obvious that 



d a w d 2 w 



dt a 

Bw = 

W(X, 0) : 



dx 2 



L-w<0 



(A.7) 



Based on (|A.7[) . next we try to prove that w = in f2y. Suppose that w obtains its 
maximum value at (x,t), if (x,t) £ Tt and the boundary conditions are Dirichlet's, 
w = holds obviously; if (x,t) 6 and u(x,t) is strictly bigger than u(x,t) for 
any (a;,t) £ f^, and the boundary conditions are Neummann's, it can be shown 
that w = holds by the ideas in Remark 14. II and the following proof. 

Now assume (x,t) 6 fix and w(x,t) > 0, then there exists t* (< i) such that 
w(x,t) < 7}w(x,i) for any t £ (0,£*) and x £ Q. Introduce the function u of t such 
that it satisfies 



du{t) 1 w(x, t) 
dt ~ 2 t* 



in (0,t*), 



and u(i) = 0, t e [i*,T], denote 



ui(a:, i) = t) + u(t) in O x [0,T]. 



Then the maximum of w is still at and w satisfies the following inequality 



d a w d 2 w T d a u(t) 
<Lw-> w 



dt a dx 2 ~ dt a 
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At (x,t), it follows that d "^t ,t] \ t =i > 0, and 
d a u{t) _ 1 



d 2 w(x,i) 

Ox 2 \x—x 



=* > 0. While 



dt° 



r(l-a) 

1 



o 



(i-r) 



t* 



< 



t* 



r a dr 



Since t < T, if T < ( 



2r(l - a) 

u>(i, t) 
2r(l - a)t* 

w(x, t) 
2r(l - a)t* 

= wpM) - a 
2r(l - a) ' 

1 \j then T Q < 1 
2LT{l-a)> ' LIiLn J 2Lr(l-a) 

M4.*) + -^ i <o. 



<It 



{t-T)- a dT 



and 



We arrive at a contradiction, so that w = holds. 



If (2lr7T^))° < r < 2(j 



: , we take T = ( 



^2LT(1-q)^ ^ - 1 — ^\2Lr(l-a)' 

O x [0, f ], we can get w = 0. For f < t < T, set t 
when 0<t<T-f, w satisfy 



- ) a . In the domain 



2LT(l-a) ' 

t - f and tu(i) = w(? + f ). So 



- L • w < 



Since T -T < ( 7 



w(a;,0) = 0. 

f )a, we have w = in Q x [0,T - f]. Then w = 



^2ir(l-a) 

£1 x [T,T]. Consequently w — in O x [0,T]. This process can be continued for any 
finite times, so we obtain w = in Qt for any T. 

Combining with the known w\ > and u>2 > 0, u> = implies wi = and 
W2 = 0. Then we get u\ = Uj and U2 = u 2 - Setting 



Ui(x,t) 



Ri, » = (1,2), 



then u(x 7 t) = (ui,u 2 ) solves (|4.4|) . 

Next we prove the uniqueness of the solution. 

Suppose that there exists another solution u'(x,t) = 
V(x,t) < u'(x,t) < U(x,t). Obviously, uj 0) < u' k < uf ] (i 



(u^u^) which satisfies 
: 1,2). On the one hand, 



since u'(x, t) can be considered as an upper solution, according to V(x, t) < u'(x, t), 

fc_>oo Si {%■, t) — Ui < u^. On the other hand, 



we know u < u' 4 . Therefore lim 
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u'(x,t) can also be considered as a lower solution, according to u'(x,t) < U(x,t), 
we know limfc_ ) . 00 u[ k \x,t) — ui > u[. Consequently u[ — Ui. 
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